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We present a novel method, based on a multiscale approach, for detecting anisotropy signatures in the arrival 
direction distribution of the highest energy cosmic rays. This method is catalog independent, i.e. it does not 
depend on the choice of a particular catalog of candidate sources, and it provides a good discrimination power 
even in presence of contaminating isotropic background. We present applications to simulated data sets of events 
corresponding to plausible scenarios for events detected, in the last decades, by world-wide surface detector-based 
observatories for charged particles. 



1. Introduction 

In many field involving data analysis, the 
search for anisotropy has played a crucial role. 
Many estimators, namely two-point correlation 
functions 0-0], have been proposed and widely 
used to search for clustering of objects and to 
measure deviation from isotropy of angular dis- 
tributions. A non-differential variant of such an 
estimators is widely used for small data set of 
points 0-0 ■ These methods apply to angular co- 
ordinates of objects as well to distributions of ar- 
rival directions of events: in this work we will 
indifferently refer to both as arrival direction dis- 
tributions of events. 

Recently, new estimators have been introduced 
to study the anisotropy signature of sky's arrival 
direction distributions: the modified two-point 
Rayleigh0, and shape-strength method derived 
from a principle component analysis of triplets 
of events [ll|. Such a test statistics have been 
recently applied to both P. Auger data and to 
synthetic maps of events, the latter generate d by 
sampling the Veron-Cetty & Veron catalog [12J 
of nearby candidate active galactic nuclei (AGN) 
within 75 Mpc (z < 0.018), showing an higher 
discrimination power than other estimators [13| . 



Within the present work, we introduce a new 
fast and simple method for anisotropy analysis, 
which makes use of a multiscale approach and de- 
pends on one parameter only, namely the angular 
scale of the instrinsic anisotropy. The main ad- 
vantage of our estimator is the possibility to ana- 
lytically treat the results: the analytical approach 
drastically reduces computation time and makes 
available the possibility of applications to very 
large data sets of objects. We test the method on 
several simulated isotropic and anisotropic arrival 
direction distributions (mock maps) and perform 
an extensive analysis of its statistical features un- 
der both the null and the alternative hypotheses. 
However, it is worth remarking that the scope 
of applicability of our method is not limited to 
UHECR physics, and it is valid for any distribu- 
tion of angular coordinates of objects. 



2. Multiscale Autocorrelation Function 

Let <S be a region of a spherical surface and let 
Pi (4>, 9) (i = 1, 2, n) be a set of points locating 
n arrival directions on S, defining a sky. The 
sky iS is partitioned within a grid of N equal- 
area (and almost-equal shape) disjoint boxes Bk 
(k = 1,2, N) as described in Ref. Q. Let O 
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be the solid angle covered by <S, whereas each box 
Bk covers the solid angle 
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where 20 is the apex angle of a cone covering the 
same solid angle: N, and are deeply related 
quantities that define a scale. 

Let ipk(®) be the density of points in the data 
set that fall into the box the function A(Q) 
that quantifies the deviation of data from an 
isotropic distribution at the scale 0,is chosen to 
be the Kullback-Leibler divergence 15, 16| 
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where ip k (Q), generally a function of the domain 
meshing, is the expected density of points isotrop- 
ically distributed on S falling into the box Bk- 
The Kullback-Leibler divergence is an informa- 
tion theoretic measure widely used in hypothe- 
sis testing and model selection criteria, statistical 
mechanics, quantum mechanics, medical and eco- 
logical studies (see [l?} and Ref. therein). This 
measure quantifies the error in selecting the den- 
sity ip(Q) to approximate the density ijj(Q): it can 
be shown that minimizing the Kullback-Leibler 
divergence is equivalent to maximum likelihood 
estimation 17j . It is straightforward to show that 



A(Q) is minimum for an isotropic distribution of 
points, or, in general, when ip(&) ~ i/)(Q), i.e. if 
the model is correct. 

If Aiata(0) an d ^iso(@) refer, respectively, to 
the data and to an isotropic realization with the 
same number of events, we define multiscale au- 
tocorrelation function (MAF) the estimator 
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where (Aj so (0)) and er.A i50 (6) arc the sample 
mean and the sample standard deviation, respec- 
tively, estimated on several isotropic realizations 
of the data. If % denotes the null hypothesis of 
an underlying isotropic distribution for the data, 
the chance probability at the angular scale 0, 
properly penalized because of the scan on 0, is 



the probability 

P (Q) = Pr( Siso (e') > .s d ata(e)|^o,ve' e v) (3) 

obtained from the fraction of null models giving 
a multiscale autocorrelation, at any angular scale 
0' in the parameter space V, greater or equal 
than that of data at the scale 0. The null hy- 
pothesis is eventually rejected in favor of the al- 
ternative Tii = -1H0 — being -i the negation op- 
erator — at the angular scale 0, with probability 

i-p(e). 

Under the null hypothesis Tio, the estimator 
s(0) follows an half-gaussian distribution, inde- 
pendently on the value of the angular scale and 
on the number of events on S 17 1 . 



3. Dynamical Boxing 

The simplest definition of the boxing algo- 
rithm, as shortly described in Sec. [2J involves 
the fixed grid introduced in Ref. [l4[ , where each 
box only embodies the relative number of events 
falling in it. Unfortunately, such a static binning 
approach could not reveal an existing cluster. In- 
deed, the fixed grid may cut a cluster of points 
within one or more edges, causing a further loss 
of information at the angular scale under inves- 
tigation. To overcome this possible loss of infor- 
mation, we introduced the dynamical binning, a 
type of smoothing of the grid by applying it on 
the data 



17] 



The smoothing, adopted in our study, deals 
with a new counting procedure for the estimation 
of the density ^(0). However, such a density 
depends on the observatory's exposure, generally 
a function on the celestial sphere depending on 
both the latitude of the experiment and the max- 
imum zenith of detection, quantifying the effec- 
tive time-integrated detection area for the flux of 
particles from each observable sky position. The 
relative exposure u> is the dimcnsionlcss function 
corresponding to the exposure normalized to its 
maximum value. For a single full-time operating 
detector, i.e. with constant exposure in right as- 
cension, fully efficient for particles arriving with 
zenith angles smaller than # max , it still depends 
on the declination 5 [3]. The proposed method 
takes into account the effect of non-uniform ex- 
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posure on the arrival direction distribution of ob- 
jects, by weighting the angular region around a 
given direction with the local value of the expo- 
sure. Our numerical studies show that such a 
dynamical binning approach recovers the correct 
information on the amount of clustering in the 
data [13]. 

4. Interpretation of MAF 

Any catalog-independent method provides in- 
formation about the angular scale 0* where the 
significance is minimum. In the case of a sim- 
ple two-point method, such an angular scale is 
quite difficult to interpret and topologically dif- 
ferent configurations of events lead to the same 
significance. In the case of the modified two- 
point Rayleigh method, the estimation of the sig- 
nificance includes another set of parameters, in- 
dependent from the angular distribution, as de- 
scribed in Ref. [Io| : parameters are sensitive to 
the orientation of the pairs and therefore to skies 
showing preferential directions and filamentary 
structure of points. It follows that O* is the most 
significant angular size for the mix of these infor- 
mations, still linked to the pair configuration. In 
the case of the shape-strength method, the esti- 
mation of two parameters, namely the shape and 
strengh, is performed: both can be interpreted, 
respectively, in terms of size and elongation of 
the triangles defined by a triplet of points. It fol- 
lows that all information is recovered from the 
configuration of triplets. 

In the specific case of MAF, the angular scale 
0*, where the significance is minimum, turns to 
be the significative clustering scale: it is the scale 
at which occurs a greater accumulation of points 
respect to that one occuring by chance, with no 
regard for a particular configuration of points, 
e.g. doublets or triplets. Fig. [T] shows the chance 
probability, as estimated by MAF, for three dif- 
ferent point source mock skies contaminated by 
80% isotropic background. In each sky, 20% of 
events are normally distributed, with dispersion 
p, around a single source: for the results shown 
in Fig.[TJ we used three representative values for 
the dispersion, namely p = 4°, p = 10° and 
p = 25°. Although the high isotropic background, 



the MAF method, by means of the dynamical bin- 
ning, is able to recover the correct dispersion: in 
all cases, the most significant angular scale for 
clustering, indicated by arrows, recovers the dis- 
persion of the corresponding sky. 
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Figure 1. MAF: chance probability versus the an- 
gular scale for three mock skies of 60 events. In 
each sky, 20% of events are normally distributed, 
with dispersion p, around a single source, and 
80% of events are isotropically distributed. 



5. Statistical analysis of MAF 

In this section, we investigate the statistical 
features of MAF by inspecting its behavior un- 
der both the null or the alternative hypothesis. 
In particular, we estimate the significance a (or 
Type I error), i.e. the probability to wrongly re- 
ject the null hypothesis when it is actually true, 
and the power 1 — (where /3 is known as Type 
II error), i.e. the probability to accept the alter- 
native hypothesis when it is in fact true. 

Null hypothesis. We generate isotropic maps of 
10 5 skies, by varying the number of events from 
20 to 500: for each sky in each map, we estimate 
the MAF for several values of the angular scale 9. 
Hence, we choose the value of G = 0* where the 
chance probability is minimum, as the most sig- 
nificant clustering scale: £>(©*) = argmine J>(0), 
properly penalized because of the scan on the 
parameter O, according to the definition in Eq. 
([3]) . Indipendently on the number of events in the 
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mock map, wc find an excellent flat distribution 
of probabilities p(Q*), as expected for analyses 
under the null hypothesis Ho, i.e. the MAF is 
not biased against Ti$. 

Because of the definitions in Eq. ([TJ, ([2]) and 
for the central limit theorem, an half-normal dis- 
tribution is expected for the estimator s(0). We 
find an excellent agreement between the distribu- 
tion for Montecarlo realizations and the expected 
one. It follows that the (unpenalized) probability 
to obtain by chance a value of the MAF, greater 

or equal than a given value So, is just 1 — erf > 
being erf the standard error function, indepen- 
dently on the angular scale 17j . Although this 



nice feature of the MAF estimator, generally the 
distribution of s max = max{,s(0)} is of interest 
for applications, because of the required penal- 
ization due to the scan over the parameter G. 
Hence, it is important to identify the distribution 
of the penalized probability p(0), if any. Intrigu- 
ingly, our numerical studies show that such a dis- 
tribution exists and it corresponds to one of the 
limiting densities in the extreme value theory (see 
[r^ and Ref. therein). In particular, the proba- 
bility density of maxima is known as the Gumbel 
distribution OH)]: 



g(x) = — exp 
a 



exp 



/' 



(4) 



where /i and a are the location and shape param- 
eters, respectively. 

In Fig. [2] are shown the probability densities of 
Smax for n = 40, 60, 80, 100 and 500 events: inde- 
pendently on n, each density is in excellent agree- 
ment with the Gumbel distribution of extreme 
values, for the parameters = 1.743 ± 0.002 and 
(7 = 0.470 ±0.002. Such a parameters correspond 
to the mean and to the standard deviation of the 
distribution, fx w 2.00 and a « 0.59, respectively 
(see [Hi and Ref. therein). It follows that the 
probability to obtain a maximum value of s(0), 
at any angular scale 0, greater or equal than a 
given value max{s(0)} is 



P (s max ) = 1 - exp 



exp 



/' 



providing an analytical expression for the penal- 
ized probability defined in Eq. ([3]). 




max[s(6)] 



Figure 2. MAF: distributions of max{s(0)} for 
n = 40,60,80,100 and 500 events. Solid line 
correspond to the least-square fit of the Gumbel 
density with parameters fi — 1.743 ± 0.002 and 
a = 0.470 ± 0.002 (x 2 / ndf = I- 1 x !0 -5 )- 



Alternative hypothesis. In order to investigate 
the behavior of MAF under the alternative hy- 
pothesis of an underlying anisotropic distribution 
of objects, we generate anisotropic maps of 10 4 
skies, by varying the number of events from 20 
to 100. In general, the anisotropy of a sky de- 
pends on several factors: for instance, in the case 
of cosmic rays, it depends on the distribution of 
sources, on magnetic fields and on propagation 



effects as energy loss or the GZK cutoff [2lj, [22 1 
(and Ref. therein). Thus, a more complicated 
approach is required for the Montecarlo realiza- 
tion of the maps. In order to estimate the power 
of MAF, we build reasonable anisotropic maps re- 
flecting in part the real- world scenario, keeping in 
mind that our purpose is to build an anisotropic 
set of events for statistical analysis and not to 
generate events mimicking real data sets with the 
best available approximation. We proceed as fol- 
lows: 



Catalog of candidate sources. Although sev- 
eral models for production mechanisms of 
UHECR are available 0, [H| (and Ref. 
therein), [23W30I ] . it is generally accepted 
that the candidate sources are extragalac- 
tic and trace the distribution of luminous 
matter on large scales 3l|. In particu- 



lar, it has been shown that correlation with 
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possible high redshift sources is unlikely 



321. whereas compact sources are favored 
33l . |34| : the recent result reported by the 



P. Auger Collaboration experimentally sup- 
ports the latter claim, showing an high cor- 
relation between the observed data and the 
distribution of nearby active galactic nuclei 



(AGN) [35l |36[. For these reasons, we use 
the Palermo Swift-BAT hard X-ray cata- 
logue of AGN with known redshift within 
200 Mpc (z < 0.047) as the reference 
catalog of candidate sources providing the 
most complete and uniform all-sky hard X- 
ray survey up to date. 

2. Source effects. Events, from each source 
in the reference catalog, are generated by 
weighting for the source flux and for the ex- 
pected geometrical flux attenuation. Hence, 
the number of events coming from a source 
is proportional to its flux and to the fac- 
tor z~ 2 : because of the small scales and 
the high energy of cosmic rays involved in 
anisotropy studies (E > 4.0 x 10 19 EeV), we 
assume a flat universe with zero cosmologi- 
cal constant (O = 1, A = 0) and nonevolv- 
ing source. Indeed, we naively take into ac- 
count the possible deflections of the parti- 
cles, due to the random component of the 
magnetic field, by producing arrival direc- 
tions gaussianly-distributcd with dispersion 
p around the source. It is worth remarking 
that such a dispersion is strictly related to 
both the injection energy and the mass of 
the particle, as well as other physical quan- 
tities dU. 

3. Background. We take into account the pos- 
sibility for a contaminating isotropic back- 
ground of the anisotropy signal, by gener- 
ating a number of events isotropically dis- 
tributed, corresponding to a fraction /j so of 
the whole data set. 

4. Detection effects. As previously explained, 
the number of events detected by a single 
fully efficient and full-time operating sur- 
face detector, depends on its own relative 
exposure. In order to take into account such 



a detection effect, we generate the events ac- 
cording to the relative exposure of the single 
detector. 

However, for a more realistic distribution of 
events, several more constraints, in general based 
on further assumptions or debated models, are 
required: the mass of the particle, the injection 
spectrum of the source, the intervening magnetic 
field, to cite some of the most important. In our 
study of the MAF discrimination power, we fix 
p = 3° , as the mean angular deviation of UHECR 
in the galactic and extra-galactic magnetic field, 
and a background fraction /j so = 0.3. 

In order to produce a likely map of UHECR, 
we choose to generate events distributed in the 
whole sky, according to the number of events col- 
lected by surface detectors in the last decades. 
In particular, we consider events with energy 
E > 4.0 x 10 19 EeV and error on the arrival di- 
rection smaller than 5°, as detected at the Sidney 
University Giant Airshower Recorder (SUGAR), 
Akeno Giant Air Shower Array (AGASA), Hav- 
erah Park, Volcano Ranch, Yakutsk and P. Auger 
Observatory (we refer to [l7| for a comprehen- 
sive list of references for the data released by 
each experiment). However, the fluxes of par- 
ticles as measured by those experiments do not 
agree each other in the absolute fluxes, and a 
rescaling is needed [Hj]. By assuming that the 
spectrum reported by the HiRes Collaboration 
[39( corresponds to the correct energy scale, the 
rescaling, based on some specific characteristics 
of the UHECR spectrum, fixes the energy shift 
factors A for the other experiments 0, [38| . Po- 
sitions, maximum zenith angles # ma x, exposures 
and energy shift factors are reported in Tab. I, 
for each experiment, as well as the number of de- 
tected events with rescaled energy E 1 > 4.0 x 10 19 
EeV (£" = XE). In Fig. [3] is shown the relative 
geometrical exposure of each single detector listed 
in Tab. I, as well as the joint exposure of all ex- 
periments. For reference, in Fig. [4^, is shown 
the all-sky data set of 102 detected events with 
rescaled energy E' , superimposed on the distri- 
bution of AGN within 200 Mpc from the refer- 
ence catalog, whereas in Fig. 0p is shown the 
mock map of simulated events according to phys- 
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ical constraints previously described. 
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Table 1 

Surface detectors: positions, maximum zenith an- 
gles # m axj exposures, energy shift factors and 
number of detected events with rescaled energy 
E 1 > 4.0 x 10 19 EeV. 
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Figure 3. Relative geometrical exposure of each 
single detector listed in Tab. I (lines and points), 
and the joint exposure of all experiments (solid 
line) . 



In Fig. [5] we show the power 1 — j3 vs. the 
number of anisotropically distributed events, gen- 
erated as described above. A sky is labelled as 
anisotropic if, for a fixed value of the significance 
a, the penalized chance probability as defined in 
Eq. ([3]) is lesser or equal than a, i.e. if the con- 
dition p(Q*) = argminep(6) < a holds for some 
angular scale 9*. In Fig. [5]is shown the power for 
two values of the significance threshold, namely 
a = 0.1% and a = 1%, estimated through the 
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Figure 5. MAF power vs. the number of events 
sampled from anisotropic mock maps generated 
as described in the text, for values of the signifi- 
cance corresponding to a = 0.1% and a = 1%. 



analytical approach. For applications, a power 
of 90% is generally required: under this thresh- 
old the method could miss to detect an exist- 
ing anisotropy signal. In the case of the MAF, 
and for the considered anisotropic mock map, 
the power increases with the number of events n 
and it is able to detect the anisotropic signal for 
n > 60, with significance a = 1%. However, by 
decreasing the significance for the statistical test, 
the power requires a greater number of events to 
reach the 90% threshold, as expected: our test 
clearly shows that the MAF provides an excel- 
lent discrimination power for n > 80. Indeed, we 
verified the agreement between the analytical and 
the Montecarlo estimations of the discrimination 
power. 

6. Discussion and conclusion 

We introduced a new statistical test, based on a 
multiscale approach, for detecting an anisotropy 
signal in the arrival direction distribution of 
UHECR, that makes use of an information 
theoretical measure of similarity, namely the 
Kullback-Lciblcr divergence, and of the extreme 
value theory Within the present work we showed 
that our procedure is suitable for the analysis 
of both small and large data sets of events, by 
applying it on several Montecarlo realizations of 
isotropic and anisotropic synthetic data sets, cor- 
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(a) UHECR events and candidate sources. (b) Mock map. 



Figure 4. a) All-sky data set of 102 detected events with rescaled energy E' > 40 EeV (see the text 
for further information) superimposed on the distribution of AGN with known redshift (z < 0.047) from 
the Palermo SWIFT-BAT hard X-ray catalogue; b) corresponding sources (squares) and smoothed mock 
map, generated for the statistical analysis (colour indicates the number of events). Equatorial coordinates 
arc shown. 



responding to plausible scenarios in the physics of 
highest energy cosmic rays. In fact, for small data 
sets as well as for larger ones, the method is able 
to recover the information about the most signif- 
icant angular scale of clustering in the data, even 
in presence of strong isotropic contamination. 

The advantages of our approach over other 
methods are multiples. First, the method allows 
an analytical description of quantities involved in 
the estimation of the amount of anisotropy signal 
in the data, avoiding thousands of Montecarlo re- 
alizations needed for the penalizing procedure of 
results and drastically reducing the computation 
time. Second, the method allows the detection 
of a physical observable, namely the clustering 
scale, in the case of a point source. In the case 
of multiple sources, the information is about the 
most significant clustering scale(s), according to 
source distribution. Third, the method is unbi- 
ased against the null hypothesis and it provides 
an high discrimination power even in presence of 
strong contaminating isotropic background, for 
both small and large data sets. Although in this 
work we referred to UHECR physics for our ap- 
plications, it is worth remarking that the method 



is suitable for the detection of the anisotropy sig- 
nal in each data set involving a distribution of 
angular coordinates on the sphere, and it can be 
adapted to non-spherical spaces by properly re- 
defining the dynamical boxing algorithm. 
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